### Beginning ####
pacman::p_load(lme4, arm, haven, stargazer, dplyr)

if(!require(broom.mixed)){
  install.packages("broom.mixed")
} else(library(broom.mixed))

get_fixed <- function (element) {
  broom.mixed::tidy(element) %>%
    dplyr::filter(effect == "fixed")
}

get_vary_any <- function (element, regex, groupv) {
 est <- broom.mixed::tidy(element) %>% 
    dplyr::filter(grepl(regex, term)) %>% 
    dplyr::filter(group == groupv) %>%
    dplyr::select(estimate) %>% 
   dplyr::mutate(estimate = round(estimate, 2))
 print(est)
  if(nrow(est) ==0){
    NA}else{
      est
    }
}

get_groups <- function(model) {
  sapply(ranef(model),nrow)
}

# +++++++++++++++++++
# Study 1 Models ####
# +++++++++++++++++++
load("./data_clean/round_1_multilevel.rdata")

rand_id <- unlist(sapply(all_models, get_vary_any,
                         "^sd__\\(Intercept\\)", groupv = "id"))
rand_narr <- unlist(sapply(all_models, get_vary_any,
                           "^sd.*Int", "narrative_num"))
vary_true <- unlist(sapply(all_models, get_vary_any,
                           "sd__StoryTrue$", "id"))
vary_AOT<- unlist(sapply(all_models, get_vary_any, "sd__AOT$", "narrative_num"))
vary_CRT <- unlist(sapply(all_models, get_vary_any, "sd__CRT$", "narrative_num"))
vary_rus <- unlist(sapply(all_models, get_vary_any,
                          "sd__RussianOrientation$", "narrative_num"))
vary_rus <- unlist(sapply(all_models, get_vary_any,
                          "sd__RussianOrientation$", "narrative_num"))
vary_anti <- unlist(sapply(all_models, get_vary_any,
                          "sd__RussianOrientation$", "narrative_num"))
vary_crt_rus <- unlist(sapply(all_models, get_vary_any,
                          "sd__CRT:RussianOrientation$", "narrative_num"))
vary_aot_rus <- unlist(sapply(all_models, get_vary_any,
                          "sd__AOT:RussianOrientation$", "narrative_num"))
vary_crt_anti <- unlist(sapply(all_models, get_vary_any,
                              "sd__CRT:AntiEuropeOrientation$", "narrative_num"))
vary_aot_anti <- unlist(sapply(all_models, get_vary_any,
                              "sd__AOT:AntiEuropeOrientation$", "narrative_num"))
num_groups <- lapply(all_models, get_groups)
id_num <- unlist(lapply(num_groups, function(x) x[1]))
narr_num <- unlist(lapply(num_groups, function(x) x[2]))

for(i in 1:length(all_models)) {
  class(all_models[[i]]) <- "lmerMod"
}

# convergence string
converg <- rep(NA, length(all_models))
for(i in 1:length(converg)) {
  converg[i] <-  ifelse(
    is.null(all_models[[i]]@optinfo$conv$lme4$messages), # I assume convergence succeeded if NULL
    "Success",
    all_models[[i]]@optinfo$conv$lme4$messages
  )
}
for(i in 1:length(converg)) {
  converg[i] <- ifelse(stringr::str_detect(converg[i], "fail"), "Fail", "Success") # if finds word fail print Fail
}

mod_stats <- rbind(
  c("Demographic Controls", rep(c("No", "No", "No", "Yes"), 4)),
  c("\\hline \\\\[-1.8ex]", rep("", 16)),
  c("Num Respondents", unname(id_num)),
  c("sd(Respondent)", unname(rand_id)),
  c("sd(Respondent, True)", unname(vary_true)),
  c("Num Stories", unname(narr_num)),
  c("sd(Narrative)", unname(rand_narr)),
  c("sd(Narrative, CRT)", unname(vary_CRT)),
  c("sd(Narrative, AOT)",unname(vary_AOT)),
  c("sd(Narrative, Russian Orientation)", unname(vary_rus)),
  c("sd(Narrative, anti-Europe Orientation)", unname(vary_anti)),
  c("sd(Narrative, Russian Orientation, CRT)", unname(vary_crt_rus)),
  c("sd(Narrative, Russian Orientation, AOT)", unname(vary_aot_rus)),
  c("sd(Narrative, anti-Europe, AOT)", unname(vary_crt_anti)),
  c("sd(Narrative, anti-Europe, CRT)", unname(vary_aot_anti)),
  c("Convergence", converg))

stargazer(all_models, type="latex", 
          title = "Study 1 Multilevel Models", 
          dep.var.labels = "Belief",
          # ^  Notice M2 is called twice. I'm going somewhere with this.
          # Below: manually supply tidied coefficients and standard errors
          # coef =  ,
          # se = list(tidyM1$std.error, tidyM2$std.error, tidyM3$std.error),
          # Omit model statistics by default...
          #omit.table.layout = "s",
          # ...but supply your own that you created (with random effects)
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          omit = c("Age", "Female", "Education", "Income"), 
          column.sep.width = "-10pt",
          intercept.bottom = FALSE,
          font.size = "tiny",
          digits.extra = 1,
          table.layout = "=ldc-ta-s-n",
          digits = 2,
          omit.stat = c("aic", "bic"),
          out = "./tables/lmer1_1.tex")

# mod.names <- vector("character", length(all_models))
# for(i in 1:length(all_models)) {
#   mod.names[i] <- paste0("(", i, ")")
# }
# stargazer(all_models, 
#           title = "Study 1 Multilevel Models", 
#           dep.var.labels = "Belief", 
#           column.labels = c(mod.names),
#           add.lines=list(
#             c("Demographic Controls", rep(c("No", "No", "No", "Yes"), 4)
#             ),
#             c("Converged", converg)
#           ),
#           type = "latex",
#           omit = c("Age", "Female", "Education", "Income"), 
#           column.sep.width = "-10pt",
#           intercept.bottom = FALSE,
#           font.size = "tiny",
#           digits.extra = 1,
#           table.layout = "=ldc-ta-s-n",
#           digits = 2, 
#           out = "./tables/lmer1.tex")

rm(all_models)

# +++++++++++++++++++
# Study 2 Models ####
# +++++++++++++++++++
load("./data_clean/round_2_multilevel.rdata")

for(i in 1:length(all_models)) {
  class(all_models[[i]]) <- "lmerMod"
}

rand_id <- unlist(sapply(all_models, get_vary_any,
                         "^sd__\\(Intercept\\)", groupv = "id"))
rand_narr <- unlist(sapply(all_models, get_vary_any,
                           "^sd.*Int", "narrative_num"))
vary_true <- unlist(sapply(all_models, get_vary_any,
                           "sd__StoryTrue$", "id"))
vary_AOT<- unlist(sapply(all_models, get_vary_any, "sd__AOT$", "narrative_num"))
vary_CRT <- unlist(sapply(all_models, get_vary_any, "sd__CRT$", "narrative_num"))
vary_rus <- unlist(sapply(all_models, get_vary_any,
                          "sd__RussianOrientation$", "narrative_num"))
vary_rus <- unlist(sapply(all_models, get_vary_any,
                          "sd__RussianOrientation$", "narrative_num"))
vary_anti <- unlist(sapply(all_models, get_vary_any,
                           "sd__RussianOrientation$", "narrative_num"))
vary_crt_rus <- unlist(sapply(all_models, get_vary_any,
                              "sd__CRT:RussianOrientation$", "narrative_num"))
vary_aot_rus <- unlist(sapply(all_models, get_vary_any,
                              "sd__AOT:RussianOrientation$", "narrative_num"))
vary_crt_anti <- unlist(sapply(all_models, get_vary_any,
                               "sd__CRT:AntiEuropeOrientation$", "narrative_num"))
vary_aot_anti <- unlist(sapply(all_models, get_vary_any,
                               "sd__AOT:AntiEuropeOrientation$", "narrative_num"))
num_groups <- lapply(all_models, get_groups)
id_num <- unlist(lapply(num_groups, function(x) x[1]))
narr_num <- unlist(lapply(num_groups, function(x) x[2]))

# convergence string
converg <- rep(NA, length(all_models))
for(i in 1:length(converg)) {
  converg[i] <-  ifelse(
    is.null(all_models[[i]]@optinfo$conv$lme4$messages), # I assume convergence succeeded if NULL
    "Success",
    all_models[[i]]@optinfo$conv$lme4$messages
  )
}
for(i in 1:length(converg)) {
  converg[i] <- ifelse(stringr::str_detect(converg[i], "fail"), "Fail", "Success") # if finds word fail print Fail
}

mod_stats <- rbind(
  c("Demographic Controls", rep(c("No", "No", "No", "Yes"), 4)),
  c("Strata Fixed Effects", rep(c("Yes", "Yes", "Yes", "Yes"), 4)),
  c("\\hline \\\\[-1.8ex]", rep("", length(id_num))),
  c("Num Respondents", unname(id_num)),
  c("sd(Respondent)", unname(rand_id)),
  c("sd(Respondent, True)", unname(vary_true)),
  c("Num Stories", unname(narr_num)),
  c("sd(Narrative)", unname(rand_narr)),
  c("sd(Narrative, CRT)", unname(vary_CRT)),
  c("sd(Narrative, AOT)",unname(vary_AOT)),
  c("sd(Narrative, Russian Orientation)", unname(vary_rus)),
  c("sd(Narrative, anti-Europe Orientation)", unname(vary_anti)),
  c("sd(Narrative, Russian Orientation, CRT)", unname(vary_crt_rus)),
  c("sd(Narrative, Russian Orientation, AOT)", unname(vary_aot_rus)),
  c("sd(Narrative, anti-Europe, AOT)", unname(vary_crt_anti)),
  c("sd(Narrative, anti-Europe, CRT)", unname(vary_aot_anti)),
  c("Convergence", converg))

mod.names <- vector("character", length(all_models))
for(i in 1:length(all_models)) {
mod.names[i] <- paste0("(", i, ")")
}

stargazer(all_models, 
          title = "Study 2 Multilevel Models", 
          dep.var.labels = "Belief", 
          column.labels = mod.names, 
          add.lines = lapply(1:nrow(mod_stats), function(i) unlist(mod_stats[i, ])),
          type = "latex",
          omit = c("Age", "Female", "Education", "Income", "OblastStrata"), 
          column.sep.width = "-12pt",
          intercept.bottom = FALSE,
          font.size = "tiny",
          digits.extra = 1,
          table.layout = "=ldc-ta-s-n",
          digits = 2,
          omit.stat = c("aic", "bic"),
          out = "tables/lmer2_2.tex")

rm(all_models)

